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Abstract 

A novel explicit and implicit Kinetic Streamlined-Upwind Petrov Galerkin (KSUPG) scheme is presented for hyper¬ 
bolic equations such as Burgers equation and compressible Euler equations. The proposed scheme performs better 
than the original SUPG stabilized method in multi-dimensions. To demonstrate the numerical accuracy of the scheme, 
various numerical experiments have been carried out for ID and 2D Burgers equation as well as for ID and 2D Euler 
equations using Q4 and T3 elements. Eurthermore, spectral stability analysis is done for the explicit 2D formulation. 
Einally, a comparison is made between explicit and implicit versions of the KSUPG scheme. 

Keywords: Kinetic Streamlined-Upwind Petrov Galerkin scheme. Hyperbolic partial differential equations. Burgers 
equation, Euler equations. 


1. Introduction 


Einite element method is one of the popular numerical methods for solving the partial differential equations nu¬ 
merically on digital computers. The standard hnite element method is well-suited to solve elliptic partial differential 
equations efficiently EElll. But this method produces oscillations for hyperbolic partial differential equations; for 
example, governing equations of convection dominated flows require additional stabilization for the standard hnite 
element based discretization. Eor such Hows, many stabilized Hnite element methods are available in the literature 
like Streamline-Upwind Petrov Galerkin (SUPG), Discontinuous-Galerkin method, Taylor Galerkin method, Galerkin 
Least-Squares method etc. The detailed discussion about these methods are given in Els]0121 El 121. Among them, 
SUPG method is one of the popular stabilized Hnite element methods used to solve high speed compressible flows 
governed by Euler equations nniini. This method introduces diffusion along the streamline direction of the flow 
which makes it stable. Apart from diffusion requirement along the streamline, SUPG method needs additional diffu¬ 
sion across the high gradient regions especially for multidimensional case. This diffusion can be controlled by using 
a shock capturing parameter which senses the shock region and adds the diffusion appropriately. Many nonlinear 
discontinuity capturing terms are available in the literature lfT3lfT4lfT5l l. 

In the Hnite volume methods, kinetic schemes (also known as Boltzmann scheme) are interesting alternatives for 
the popular Riemann solvers. Development of these schemes are based on the fact that one can recover the Euler 
equations by applying a suitable moment method strategy to the Boltzmann equation. The Boltzmann equation is 
given by 


dl 

dt 


+ v.V/ = 



( 1 ) 


where / and v are velocity distribution function and molecular velocity respectively. The right hand side is the collision 
term and left hand side consists of an unsteady term and a convection term. The well-known BGK model simplifies 
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the collision term and converts the otherwise integro-differential equation to a partial differential equation with a 
relaxation source term M- Using an operator splitting strategy by which the solution of the Boltzmann equation is 
split into a convection step and a collision step and further employing an instantaneous relaxation to equilibrium in 
the collision step leads to a simplification which is often used in Boltzmann schemes. The moments of the resulting 
Boltzmann equation then leads to Euler equations of gas dynamics, with the equilibrium distribution /“ being a 
Maxwellian. The Euler equations can then be written in the following moment form. 

|vpJg+v.V/ = oj,/ = /"^ (2) 

Here < . > is an appropriate moment and 'T is the moment function vector. One can also obtain the Burgers equation 
using same procedure stated above by dehning an appropriate equilibrium distribution function. The advantage of 
this procedure is, instead of dealing with a nonlinear hyperbolic system of equations (Euler equations) we are dealing 
with a linear scalar equation (Boltzmann equation without collision term). There are many kinetic schemes available 
in the literature like Beam scheme of Sanders & Prendergast 03, the method of Rietz El, the Equilibrium Elux 
Method of Pullin IfTSl . Kinetic Elux Vector Splitting (KEYS) method of Deshpande Il20ll24l . the compactly supported 
distribution based methods of Kaniel m and Perthame 1221 . the Peculiar Velocity based Upwind (PVU) method of 
Raghurama Rao & Deshpande ll25l and the BGK scheme of Prendergast & Xu ll2^ . These methods were developed 
in the framework of finite difference or finite volume methods. The application of finite element methods in the 
framework of kinetic schemes is of currently ongoing interest, with some of the works in this category being due to 
Deshpande & Pironneau ll26l and Deshpande Il27l . Yu & Dai lOOll . Khobalatte & Leyland Il28ll . Tang & Warnecke 
El, Liu ans Xu |[33, Ren et al. El, Gassner Il29l . 

In this paper an attempt has been made to take advantage of the strategy of kinetic schemes for developing an 
efhcient SUPG scheme in the framework of Boltzmann schemes. Along with this novel scheme (KSUPG), we have 
also developed a simple shock capturing parameter which senses the jump inside the element for 2D Euler equations. 
However, unlike the traditional SUPG method, the shock capturing parameter is needed only for 2D Euler equations 
(not in one dimension) and not even for 2D Burgers equation. Constructing the stabilization parameter r (which is 
the intrinsic time scale) in mutidimensional SUPG framework is not a trivial task. Many methods are available in 
the literature El El- But, in the proposed KSUPG scheme r is defined for both scalar and vector equations simply 
from the linear scalar formulation. The efficiency of the new scheme is demonstrated by solving various test cases. 
This paper is arranged as follows. In section 2 governing equations for high speed flows (Euler equations) and scalar 
Burgers equation are given. Section 3 and 4 give the ID and 2D explicit KSUPG weak formulation for both Burgers 
equation and Euler equations. In section 5, a simple shock capturing parameter is introduced. Section 6 explains 
the spectral stability analysis for explicit KSUPG scheme. An implicit KSUPG formulation for ID and 2D Euler 
equations is given in section 7 followed by section 8 where various numerical test cases for both explicit and Implicit 
formulation are solved. Before ending section 8, comparison for explicit and implicit KSUPG schemes is made based 
on the number of iterations required to bring down the residue below the specified tolerance limit, the computational 
cost and the sparsity pattern of the coefficient matrix of global system of equations. 


2. Governing Equations 


The governing equations are for the inviscid compressible flows, given by Euler equations as 


dt dxi 


€ Q X [0, T] 


(3) 


where U = [p, pu\, pu2, pu^, pE]^ is the conserved variable vector and G/ = [pu,, dnp + pu\Ui, dap + pu2Ui,6aP + 
pu^Uj, pui + pujE]^ is the inviscid flux vector, p, u\,U2, U2,E, p are density, velocity components in x,y and z direc¬ 
tions, total energy and pressure respectively and 6ij is a kronecker delta. Total energy is given by 


E = 


piy- 


-h -(u, + ui 

1 ) 2 ^ ‘ 2 


ui) 


2 



Note that A, = ^ is an inviscid flux Jacobian matrix for the domain Q e M® (where D is the spatial dimension) 
with boundary F = F^ U Fat and final time is given by F e As the eigenvalues of A, are real and eigenvectors 
are linearly independent, the system of equations is hyperbolic. Beyond being hyperbolic, these equations are non¬ 
linear and produce shock waves, expansion waves and contact discontinuities which need to be resolved in numerical 
simulations. 

We also consider a scalar hyperbolic conservation law as 


^ dgiju) 
dt dxi 


e Ox[0, T] 


(4) 


where u is the conserved variable. The fluxes gi(u) can be linear or nonlinear. One example is the inviscid Burgers 
equation in the the fluxes are nonlinear and produce shock waves and expansion waves. 


3. One Dimensional Explicit KSUPG Weak Formulation 

The standard Galerkin finite element approximation for molecular velocity distribution function is 

(5) 

Vi 


where the domain is divided into riei elements. 


tlel 

O = O,' and O,- n Oy = 0, Vi j 

Vi 


( 6 ) 


Defining the appropriate test and trial functions spaces as = {N^ e “TF' 

and = 0 on F/j) and = {/^ e FF' and /* = /^ on Fo) where F/j is the Dirichlet boundary, the weak formulation 
is written as, find e such that V A* e 



A" H- TV- 


c/Af 

dx 


dQ.i - 0 


(7) 


where t = h/(2|v|). The global system of equations are obtained as 




h dN^\ 
2|v| dx ) 


dQ = 0 


or 




( 8 ) 


where basis functions A^ 6 C*'(0). It is important to note that, the test function is enriched with additional term which 
is multiplied only with the convection term. That gives a required diffusion term. 

In matrix form. 


+ Cvf + sign(v) v/ H- /a, = 0 
at 2 


(9) 
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where Mass matrix M, Convection matrix C, Diffuion matrix D and Neumann boundary condition are given by 



All the integrals are evaluated with full Gauss-Quadrature integration. Taking moments with the suitable moment 
function vector T* 

+ Cvf + ^Dsign(y) vf + 

or +C <'V,vf > +-D <'¥, sign(y) y/ > -k /^ > = 0 (10) 

dt 2 

Equation ( [T0| i is the semi-discrete weak formulation. 

3.1. One Dimensional Burgers Equation 
The ID Burgers equation is given by 


du dg{u) 
dt ^ dx 


e Gx[0, T] 


( 11 ) 


For the sake of convenience, we write the flux g{u) = as g{u) = cu with c - |m. In case of one dimensional 
Burgers equation 'T = 1 and Maxwellian distribution function to recover the Burgers equation as a moment from the 
Boltzmann equation is given by 


/ 


M 





where c is constant to recover the linear convection equation and is a function of u, i.e., c = m/2 to recover the inviscid 
Burgers equation. For one dimensional problem D - Fet us now evaluate the terms for the case of one dimensional 
Burgers equation. 



^OO 

I dv - u 

/ —OO 

(12) 

<'T,y/" > = 

I vf^ dv — cu 

/ —OO 

(13) 


^CO 

< T*, sign(y) vf^ > - I sign(y) vf^ dv 


\J —CO 

^0 /^oo 


vf^dv+ vf^dv 


%J —CO %Jo 


c, ^ “ -s- 

= CM erf(i) H — —e 
sJtiP 

(14) 

r 5 < T', f'* > r du 

<'V,fN>^ - 

JYn JYn 

(15) 
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where s - u ■v^/2 and /3 - I- Note that, since no energy equation is involved (no pressure and temperature terms) so, 
/? is just a constant value calculated from the standard Maxwellian distribution function. Moments of last expression 
lead to the Neumann boundary condition in macroscopic variable u. 

Substituting these values in equation ( [TOl l, we get 

du hi M 

M -h CcM H—D cMerf(i)H- —e + mat = 0 (16) 

dt 2 \ j 

where is the Neumann boundary condition in macroscopic variable. 


3.1.1. Temporal Discretization: 

In this work finite difference approach is adopted for temporal discretization using 6 method as 


M- 


At 


(1 -0) CcV + -D 


c"u’' erf(i) + 


+ 0|Cc"m" + 


erf(s) + 


sfjip 

,,n+\ 


sJnP 


+ — 0 


(17) 


Thus, 0-0 gives an explicit method and 0=1 gives an implicit method. Semi-Implicit methods can be obtained with 
0 < 0 < 1. For example, 0=1/2 gives Crank-Nicolson method. 


3.2. One Dimensional Euler Equations 
The ID Euler equations are given by 


where 


U = 


dt dx 


0 


p ) ( pu 

pu > , G = f p + pu^ 
pE j y pu + puE 


(18) 


are the solution vector and the flux vector respectively. For recovering the ID Euler equations as moments of the 
Boltzmann equation, the Maxwellian distribution function is given by 


/ 


M 




(19) 


where v is the molecular velocity, / is the internal energy variable corresponding to the non-translational degrees of 
freedom and is the average internal energy corresponding to the non-translational degrees of freedom which is given 
by 

3 — y 

lo =- —RT (20) 

2 ( 7 - 1 ) 

and 7 is the ratio of specific heats. 

For ID Euler equations, moment function vector T* is defined as 


T' = 



( 21 ) 


With the Maxwellian / = the other terms in weak formulation given by equation ([TOli can be obtained as 











( 23 ) 


X oo ^oo pu 

dl '¥vf^dv^\ p+pu^ i = G 
( pu+puE j 


/^CO /^OO 

<'F, sign(v) v/“ > = I dl I sign(v)vf^ dv 

\Jo «J —oo 


J r-»CXD W) ^CO /-*€<} 

dl j vf^dv+ j dl j vf^ dv 
0 j-oo Jo Jo 


Taking the first integral term 


f^CX! 


puA —pB 

f dl 

II 

(p + pu^)A~ — puB 

Jo 

J—oo 

{p+pE)uA--pi^^+E)B 

can be evaluated as 


/^QC 

^oo 

puA^ +pB 

( dl 

vf^dv^- 

ip + pu^)A^ + puB 

Jo 

Jo 

{p+pE)uA^+ p[-^+E)B 


where 


l±erf(s) 1 


= 


and B = 


2 2 

Here, s = m p — IjlRT. Substituting these values in equation (|24| and then simplifying we get 


< T', sign(v) vr >^< 


PM erf(i)+ 

(p + pM^) erf(s) + 

(p+p£)Merf(i)+p(^ +E)^e- 


Substituting these values in ([TOli we get 


d \ P \ 

M— f pM f + p + pu^ 

^ [ pE j [ pM + puE 


h 

+ -D{ 


putxf{s)+-^^e-^ 

(p+pM2)erf(s)+ 

(p+p£)Merf(i)+p(^+£) -^e' 


' Um 


= 0 


(24) 


(25) 


(26) 


(27) 


(28) 


3.2.1. Temporal Discretization: 

Temporal discretization is done using 6 Method with 6 = 0 as 

ij«+i - \J" i h \ 

M - — -+ CAU" + -D < T', sign(v) v/" >" 1 + m^t = 0 (29) 

In the above discretized form the flux vector G is written as AU where A is the flux Jacobian matrix given by 


A = 


0 1 0 

(^)m^ (3-y)M (y-l) 

(■y-l)M -yMii 'YE-^(y-\)u^ yu 


(30) 
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3.3. Linearization and Iterative Solver 

The global nonlinear fully discretized equation of the form K(u)u - (p can be linearized by using Picard iteration 
technique as 

= (j) 

Then, the linearized system of equations is solved using bi-conjugate gradient stabilized method. 

4. Two Dimensional Explicit KSUPG Weak Formulation 

The standard Galerkin finite element approximation for molecular velocity distribution function is 

(31) 


where the domain is divided into n^i elements. 


iiel 

G = Q,' and Q,- n Llj - 0, V/ j 

Vi 


Defining the test and trial functions spaces as 'V'’ 


{At* 


(32) 


“H' and 


- 0 on To) and - {/^ e “H' and/* = /* on To) where To is the Dirichlet boundary, the weak formula¬ 
tion is written as, find /* e such that V e 


^ r idf df df 


W* -H 


dm 


dm 


TlVi—— +TIV 2 —J— 

dx dy 


dLli - 0 


where ti = /i/(2|vi|) and T 2 = /i/(2|v2|). The global system of equations are obtained as 


r [df df df\ 
Jni dt ’’’ dx dy] 


IIA^* H- 


dm 


dm 


TlVi-— -HT1V2-— 

dx dy 


dQ^O 


(33) 


(34) 


h-xT 


(N") 


h 

^ 2 
h 


f {ff{f)dQ.% -I- r 

Jq Jq 

I(^) 

r idmf idm\ 

1 c/Gsign(vl) v2/ 

dQ sign(v2) v2f 

in\ oy I \ uy / 

X ■ 


dm 

^-\dLlvf + 


X 


h\T 


im) 


dm 

^|^0V2/ 


^ Jn \ dx ) \ dy 

r ! n\jh\^ I nAjh 

2 Jn \ «9x / \ dy } 
h r IdN^f ldN''\ 
^ 2 Jn \ 5 y j \ ( 9 y / 

dfh 

+ I = o 

Jvk dn 


(35) 


where basis functions e C°(f2). Again, the enriched terms present in the test function are multiplied only with 
convective terms which gives diffusion terms in x, y directions and cross-diffusion terms in x - y directions. 
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In matrix form, 


df h h 

+ CxV\f + CyV2f + -Z)^sign(vi) vi/ + -D^j,(sign(vi) V2/ + sign(v2) vif) 

+ ^£>ySign(v 2 ) V 2 / + /iv = 0 (36) 

where 


r (N''f(N'')dQ 
Jn 



All integrals are evaluated using full Gauss-Quadrature integration. Taking moments 

I df h h 

\ ’ ^~dt ^ ’’’ 2^xy{^^E,'^{v\) V2f + sign(v2) vi/) 

h \ 

+ -DySign(v2)v2f) + <fN>^0 

^ 'l',vi/ > +Cy < V 2 / > < T',sign(vi)vi/ > 

at 2 

+ T', sign(vi) V 2 f> + < T*, sign(v 2 ) vif>) + ^Dy < T, sign(v 2 ) V 2 / > 

+ <T',/^>=0 

Now lets evaluate these moments for 2D Burgers equation and 2D Euler equations. 


4.1. 2D Burgers Equation 

The 2D Burgers equation is given by 


du 

dt 


dg\{u) dg 2 {u) 


dx 


dy 


= 0 


is written as 


du 

dt 


dciu dc2U 


= 0 


(37) 


(38) 


dx dy 

where ci and C 2 can be functions of u for obtaining nonlinearity or can be constants for keeping them as linear. 
Maxwellian distribution function for recovering the 2D Burgers Equation as a moment of the Boltzmann equation is 
given by 


f 


^“0 


e 


{-P{V\-Cif-P(V2-C2f) 


( 39 ) 



Here, *? = 1 and value of f} is fixed as unity. 

For 2D Burgers equation one can evaluate the integrals as 


<q/ yM 


> = I j dvidv2 - u 


(40) 




> = 



vi/^ dvidv2 - ciu 


(41) 


<'¥,V2f’^>- I I V2f'^dvidV2-C2U 


(42) 




I 


d < '¥,f > 

dn 


dY\r — 


X 


du 
r« dn 


< 'I',sign(vi)vi/'^ > 


< 'F, sign(v2), V2f^ > 


<'F,sign(v 2 )vi/"> 


<'F,sign(vi)v2/"> 


/-»00 /-»CXD 

si 

\J—CO \J—CO 


sign(vi)vi/"(fvit/v2 


P 

—u 

n 


e 1 

-+cierf(si) 

n 


/-»c» /-»CXD 

si 

%J—CO *J—CO 


sign(v2) V2f^ dvidv2 


P 

—u 

71 


e 2 

-+ C2erf(s2) 

71 


/-»C» /-»CXD 

si 

%J—CO *J—CO 


sign(v2)vi/"(fvit/v2 


—uc\Qrf{s2) 

71 


/-»oo /-»oo 

\)—CO k)—< 


sign(vi)v2/ c/vi<iv2 


MC2erf(si) 


(43) 


(44) 


(45) 


(46) 


(47) 


4.2. 2D Euler Equations 


The 2D Euler equations are given by 


dV dGi dG2 

dt dx dx 


where 


1 " 1 r 

f pui V, Gi = < 

pui 


pU2 

PU\U2 
, 2 

‘, G 2 — < 

p + pu\ 

( pu2pE ] 

P + pUy 
pui + puiE 


pUiU2 

PU2 + PU 2 E 


are the solution vector and the flux vectors in x and y directions resptively. 

In case of 2D Euler equations, the Maxwellian distribution function is given as 


/ 


M 


P_ P_ Jy-l3(V\-Uif-l3(V2-U2f- ) 

lo 71 


(48) 


(49) 
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where vi and V 2 are molecular velocities in x and y directions and is defined as 


2 - y 

lo = 7 

(r- 1) 


The vector T* is defined as 




1 

Vl 

Vi 

7+^ 
^ ^ 2 


and p - ^. Integrals are evaluated as 


< 'F 


^00 ^oo ^oo 

I dl \ dvi dv 2 ^f^ 

kJO kJ—OO kJ~00 


p 

pUi 

PU2 

pE 


= u 


^oo ^oo ^oo 

< T', Vi/" > = I dl j dvi dv2'^vif’^ 

*J0 \J—OO v/—OO 


pu\ 

P + pui 
pUlU2 

pui + pu\E 


= G, 


Flux vector G\ is further decompose by using homogenity property Gi = AiU where 


Ai = 


0 


1 


0 0 

+ (3 -y)Mi -(t-1)m2 7“! 

— U\U2 U2 Ml 0 

-{yE - (y — 1)(mj + u^)u\ yE — ^(2mj(Mj + m^)) —(y — l)uiU 2 yuy 


/~*0O /~*oo /~*oo 

< 'P, V 2 /" > = I dl \ dvi dV 2 'l'y 2 f^ 

%J0 *J~00 kJ~00 


pU2 
pUiU2 
p+ pU2 
pU2 + PU 2 E 


= G, 


Similarly, flux vector G 2 is further decompose by usign homogenity property G 2 = A 2 U where 


A2 = 


0 


0 1 

-U\U2 U\ W 2 

-Mj + ^Cmj+m^) -(7-1)mi (3-y)u2 

-{yE - (y - l){u] + ul))u 2 -{y - l)uiU 2 yE - ^{2u\{u\ + u^) yu 2 


0 

0 

y- 1 


<^Jn> 


-i 


8 < T',/'’ > 


dn 


N''dYj^ 


on 


( 50 ) 


(51) 


(52) 


(53) 


(54) 


( 55 ) 
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< sign(vi)vi/'^ > = < 


PMierf(si)+p^ 

(p+PMi)erf(si)+pMi^ 

_^2 

puxU2&cf{si) + pU2^ 

(^P + 5 P(“? + “2)) «ierf(ii) + (z^P + 5 P(m? + “2)) ^ 


< sign(v 2 ) V 2 f^ > = < 




pM2erf(i2)+p^ 

pMiM2erf(s2) +pMi^ 

_^2 

(P + pM2)erf(s2) + PM2 ^ 

(^P + 5P(“? + “ 2 )) «2erf(i2) + (z^P + 5 P(m? + “ 2 )) ^ 


(56) 


(57) 


<'I',sign(vi)v2/" >= 


pM2erf(ii) 
PM2(^ +Mierf(si)j 
perf(ii)(^ +m|) 


pW2 


4( 


p/oM2erf(si) + f erf(si)(^ + m^) 


2 _ 

yS\(S 


-^eif(ii) 


2 mi 


vs 


erf(ii) 


V^) 


(58) 


Similarly, 


< 'F, sign(v2)vi/'^ >= 


pui erf(i2) 
perf(i2)(^ +M^) 

PMi(^ +“2erf(i2)) 


I PU\ 
2 



p/oMieif(i2) + f erf(i2)(^ 


2 


-52^ 2 


+ ^eif(i2) 


. 2m2 

+ —-p 
-f ^ c 


+ ul) 

-i + ^erf(^2)V^) 


(59) 


As usual temporal discretization is done by using 0 method with 0-0 for explicit KSUPG scheme. The non¬ 
linear system of equations is linearized using Picard iteration technique and is solved by using bi-conjugate gradient 
stabilized method. 


5. Shock Capturing Paramter 

In case of multidimensional KSUPG method, diffusion along streamline direction is not sufficient to suppress the 
oscillations near high gradient region. Hence additional diffusion terms with a shock capturing parameter is required 
which can sense these high gradient regions and adds additional diffusion. There are many shock capturing parameters 
available in the literature Qoiiia. In this work we present a simple gradient based shock capturing parameter as 
follows. 

We define a simple element-wise gradient based shock capturing parameter which introduces diffusion along 
high gradient directions. Figure (a) shows a typical four node quadrilateral element. As shown in figure, the 
maximum change in T' (where T* could be density, temperature or even pressure; in present work, density is used 
for all numerical test cases because density is the only primitive variable which jumps across all the three waves: 
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(a) 



Figure 1: Four node quadrilateral element in physical domain 


shocks, expansion waves and contact disctonintuies) occurs across node 1 and 3. The element based shock capturing 
parameter is then defined for node 1 and 3 as 


,ele W llV'Pll^n 

I l|T'||S'= / 

where subscripts 1 and 3 represent node numbers. For nodes 2 and 4, it is defined as 


cele 

Oi = - 

a 


Vi/ele _ \i/ele \ 
^Max 

lixi/liele 


where i = 2,4. 


(60) 


( 61 ) 


Here 1.4 < cr < 2. For most of the test cases a — 2 works fine. At element level matrix form, the shock capturing 
parameter is given by 


^ele 


(5f 0 0 

0 0 

0 0 
0 0 0 


The upper and lower bound on the value ||VT'||^‘^ is given by 


0 

0 

0 


0 < iivT'ii^‘= < 


(62) 


(63) 


It is important to note that, the addition of extra shock capturing term in the weak formulation makes the formulation 
inconsistent with the original equation. Thus, we define such that as h 0, should disappear. This condition 
is achieved by including h in the numerator, which vanishes as we refine the mesh. Similarly, one can define such 
a delta parameter for triangular elements shown in figure [T](b). The additional diffusion term along with the shock 
capturing parameter is then given by 

6(D, + D,,) (64) 

where 6 is the global matrix obtained by assembly and Dy are the diffusion matrices in x,y direction respectively. 
These diffusion matrices are defined in 2D Euler KSUPG formulation. 


6. Spectral Stability Analysis 

Stability analysis of a numerical scheme gives the acceptable value of time step At within which the scheme is 
stable. In other words, error does not grow with time. Unlike von Neumann stability analysis, spectral stability 
analysis includes the boundary points too. In the following analysis, we consider the 2D weak formulation of a linear 
equation. The global system of equation can be written as 

j/f+Af ^ 
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where ^ is amplification matrix and U' are the numerical solution at time level f + Af and t respectively. Let U 

be the exact solution, then the error e is given by 

= U-U\ Vt6K+ (66) 

Substituting this in equation (|6^ one can obtain 


= Jle' = = 




(67) 

( 68 ) 


Rearranging, we get 


For stable solution 


j+ht 




•t+bj 


j+ht 


< 1 


which gives < 1 —> ||,;?1|| < 1. We use the following relation 

\gm\ < m\\ < 1 


(69) 


where gi^) is the spectral radius of amplification matrix. 

Thus, error e‘^^' remains bounded when the maximum eigenvalue of amplification matrix ^ is less than or equal 
to unity. To find ^ matrix, we use explicit weak formulation for 2D scalar linear problem which is given by 


M- 


- m" 


At 


F-\ 


e '^1 

-+cierf(5i) 



n 


+Dx3,(c2erf(si) + cierf(s2))M" + Dy 
here c\ and C 2 are constants. Rearranging above equation, we get 


e 2 

-+ C2erf(s2) 

n 


m" U 0 




M - Atc\Cx + c^Cy + - -i/- 

2 V TT 


e 1 

-+ cierf(ii) 

TT 


+Djcyic2srf(si) + cieif(i2)) + Dy 


e ‘*2 


+ C2erf(jf2) 


which gives matrix as 


Jm-A f 

+ C2Cy + - -,1^ Idjc 

r -^2 1 

e 1 

-+cierf(si) 

1 


TT 


+D„(c2erf(si) + cierf(s2)) + Dy 


-+ C2erf(s2) 

TT 


■IM 


Using the stability condition given by equation ( [69] l, we get 

At, 


e[i-M-^[cTCy + c2Cy + 'i^r)] 


(70) 
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Figure 2; Spectral radii of amplification matrix ^ at different time steps At. 


where 


-s^ 

e 1 


-s^ 

e 2 

--Hcieif(ii) 

+ D^{c2e.rf{sx) + cierf(i2)) + Dy 

-H- C2erf(i2) 

71 

71 


The maximum eigenvalue A^ax (which is the absolute maximum of all eigenvalues ) is computed numerically using 
Rayleigh quotient over a 32 x 32 grid for 2D linear convection equation with unity wavespeeds in both directions. The 
initial condition is a cosine pulse convecting diagonally in a square domain [0,1] x [0,1]. 

Figure |^shows the spectral radii of amplification matrix at different time steps At. It shows the solution become 
unstable after At - 0.005. 


7. Implicit KSUPG Formulation 

Explicit schemes are slow because of the small values of time-steps used based on the conditional stability limits. 
CFL condition puts time step restriction on such schemes. Therefore, solving problems like vector conservation laws 
{e.g., Euler equations) using explicit schemes are computationally very expensive. An alternative is to use implicit 
schemes which are unconditionally stable without any restriction on time step. But this advantage comes with extra 
liability. In this case the coefficient matrix of the solution which needs to be inverted is quite complicated, hence 
require utmost care. 

In this section, an implicit formulation of the KSUPG scheme is derived which is computationally very efficient. 
This formulation is derived for ID as well as for 2D Euler equations. 

7.1. One Dimensional Euler Equations 

Using 6 method for temporial discretization of ID Euler equations 

Tjn+l _ TJ« / U 

M - — --H (1 - 6») CAU"+^ H- -D < T', sign(v) v/" 

-I- 6»|CAU" -H < 'P, sign(y) v/" >"| = 0 


for 0 - 1, we get 


M- 


u« 


■U" 


At 


+ I CAU"+' -H -D < 'P, sign(v) v/" >"+‘ 


= 0 
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where < 'P, sign(v) vf^ >= DU 


D = 


p erf( 5 ) 

P 

_,2 


0 

pevi{s) 

P 


0 

0 


Merf(s) + 



7.2. Two Dimensional Euler Equations 

Again, using 9 method for 2D Euler equations 

Tjn+l _ IJ« / h 

M - - -+ (1 - 0) I C,G” + C,G” + -D, < 'P, sign(vi) vi/" >" 

+ '^D,y{< '¥, sign(vi) V 2 /" >" + < 'P, sign(v 2 ) n/" >") + < 'P, sign(v 2 ) V 2 /" >"j 

+ 0|c,Gf' + CyGf 1 < 'P,sign(vi) vi/" >"+' 

+ ^D,,(< 'P,sign(vi)v 2 /" >"+' + < 'P,sign(v 2 )vi/" >"+') + < 'P, sign(v 2 ) V 2 /" >"^'j = 0 

For implicit method 9-1, which gives 

Tjn+l _ Tjn / j, 

M - - -+ (c,Gf' + C,G;^' + -D, < 'P,sign(vi)vi/" 

+ \D,y{< 'P,sign(vi)v 2 /" >"+' + < 'P,sign(v 2 ) vi/" >"^') + < 'P, sign(v 2 ) V 2 /" >"^'j = 0 


(71) 


Decomposing vectors Gi = AiU, G 2 = A 2 U, < 'P, sign(vi) vif^ >- DxU, 

< 'P, sign(v 2 ) V 2 /^ >= DyU, < 'P, sign(vi) V 2 /^ >= D^yU, and < 'P, sign(v 2 ) vi/''^ >= Dy^U, where Ai, A 2 are 
Jacobian matrixes defined previously and 


Dx 


Ml erf(si) + 

r’erf(.vi) 


_2 
e 1 


•EiP 


P 


0 

0 


0 

Mierf(ii)+^ 

0 

0 


0 

0 

Ml erf(ii) + ^ 
0 


0 

0 

0 

©11 


(72) 


0 0 0 

M2 erf(s2) + ^ 0 0 

0 M2 erf(s2) + ^ 0 

0 0 022 


M 2 erf(si) 

0 

0 


0 

0 

0 

Ml erf(yi) + 

e *1 

■v^ 

0 

erf(si)(^ +m2) 

0 

0 


0 

0 

0 

0 


©12 

Ml erf(s2) 


0 

0 

0 

erf(s2)(^ +“i) 


0 

0 

0 

0 

M2 

erf(*2) + ^ 

0 

0 

0 


0 

0 

©21 


Dy = 


M2erf(s2)+^ 

0 

perf(i2) 

P 

0 


(73) 


(74) 


( 75 ) 
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with 


©n = 


[{^P + 2^“? + M2))«ierf(^i) + {^)P + 5P(“? + “ 2 )) ^ 

P£ 


©22 - 


{(^/^ + sMm? + M 2 )) M2erf(i2) + + 5P(“? + “ 2 )) ^ 

P£ 


p ^/3M2 , ,3 


©12 =ip/oM2erf(si) + +-erf(si)|— + M 2 


^ PM2 m 
2 V TT 

UpE) 


2mi 


18 # 


-iie *1 V)r , 
— + ^erf(.0 


+ - 1) H-^erf(ii) -y^ 

/? # y 


©21 = '{p/o Mierf(i 2 ) + I erf(i 2 ) + u\ 


pu\ P 


2u2 2 

+ 


2 y 7T 

lipE) 

the final implicit equation is written as 

U«+i _ u« 

M- 


P P^fP 


-S2e *2 r, . 

+ ^erf(«) 


2 m 9 2 _ _ 

+ —(e " - 1) +-|erf(s2) V)r 
P W 


At 


+ + CyAz + - [DxDj. + (Dxy + Dyx)£>xy + DyAjjU”^' - 0 


8. Results and Discussion 

In this section various test cases are solved for inviscid Burgers as well as Euler equations to demonstrate the 
accuracy, efficiency and robustness of the proposed scheme. These codes are run on 3.10 GHz desktop machine. 


8.1. 1 -D Inviscid Burgers equation test case 
Consider the inviscid Burgers equation in 1-D 


du 

dt 


dx 



= 0 


with initial conditions representing a square wave given by 


m(x, 0) 


1 for |x| < 1/3 
-1 for 1/3 < |x| < 1 


(76) 


Note that c - m/2 for this equation. 

In this test case, the jump from -1 to 1 at x = -1/3 creates an expansion fan whereas the jump from 1 to -1 creates 
a steady shock wave. The solution at time t - 0.3 contains a sonic point (where velocity m = 0) in the expansion fan. 
Many schemes generate unphysical solutions at sonic points due to insufficient numerical diffusion. Figure shows 
the numerical solution obtained with 50 grid points using CEL number 0.3. The proposed method does not encounter 
any sonic point problem. This scheme captures the shock with just two grid points. 
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Figure 3: Inviscid Burgers Equation Test Case 


8.2. Test cases for ID Euler equations 

For ID Euler equations we shall solve the following test cases. 


8.2.1. Sod’s Shock Tube Problem: 

Sod’s shock tube problem consists of a left rarefaction, a right shock wave and a contact discontinuity which 
separates the rarefaction and shock wave. The initial conditions are given by 


p(x, 0), m(x, 0), p{x, 0) 


|l, 0,100000 If -10<x<0 
|o.l25, 0,10000 If0<x<10 


(77) 


The number of node points are 100 and CEL number is 0.15. Final time is t = 0.01. Figure shows the density, 
velocity, pressure and Mach number plots. Here, all the essential features like expansion wave, contact discontinuity 
and shock wave are captured reasonably well. 


8.2.2. Shock Tube Problem of Lax: 

The initial conditions for the shock tube problem of Lax are 


p{x, 0), m(x, 0), p{x, 0) 


fo.445, 0.698,3.528 If 0 < x < 0.5 
|o.5, 0,0.571 If0.5<x<l 


(78) 


The number of node points are 100 and CFL number is 0.1. Final time is t=0.13. Figure |^shows the density, velocity, 
pressure and internal energy plots. 


8.2.3. Strong Rarefactions Riemann Problem 

The initial conditions for the Riemann problem are 


p(x, 0), m(x, 0), p(x, 0) 


fl,-0.2,0.4 If0<x<0.5 
|l,2,0.4 If0.5<x<l 


(79) 


In this test case a near vacuum state is reached. Many popular schemes like linearized Riemann solver fails to predict 
the correct pressure and density and typically give negative pressure and density. The number of node points are 200 
and CFL number is 0.1. Final time is t=0.15. Figure [^shows the density, pressure and velocity plots. 
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Figure 4: Sod’s Shock Tube Problem 


8.3. Two dimensional Burgers equation 

Two dimensional Burgers equations is given by 


du ^ ( 2 ) du 
dt dx dy 


(80) 


The boundary conditions are: 

u(0,y) - 1 and u{\,y) = -1, 0 < y < 1 

and 

u(x, 0) = 1 — 2x, 0 < X < 1 

Exact solution is given in ll35l . 

Here, c\ - ujl and C 2 = 1. Figures |^and |^show the contour and surface plots of steady state solution on 32 x 32 
Q4 and T3 meshes respectively. The normal shock wave is captured quite accurately using such a coarse grid. Figure 
l^shows the residue plot for Q4 mesh where residue is given by 


Residue = 




8.3.1. Comparison of Standard SUPG and the Proposed KSUPG Scheme 

The comparison has been done with the standard SUPG method with shock capturing parameter used in ifTSll with 
the proposed KSUPG scheme. Figure shows the result of 2D Burgers equation over 32 x 32 Q4 element grid. 
Accuracy of the solution is more in the proposed scheme even without using shock capturing parameter. 
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Figure 5: Shock Tube Problem of Lax 


8.4. 2D Euler Test Cases 

Unlike in the case of 2D Burgers equation, additional diffusion is required for simulating 2D Euler equations and 
thus the shock capturing parameter given by equation ( [64l i is added in both explicit and implicit KSUPG formulation 
for 2D Euler equations. 


8.4.1. Oblique Shock m: 

The domain is square [0, 1] x [0, 1], the left boundary is the inlet with Mach 2 at an angle of -10" to the bot¬ 
tom boundary. Bottom boundary is the wall from where oblique shock wave is generated which makes an angle of 
29.3° with the wall. The Dirichlet boundary conditions on left and top boundaries are p - \,u\ — cos 10°, M 2 = 
- sin 10°, p — 0.179. At the wall, no-slip condition is applied, i.e., \.n - 0 where v is a velocity vector in two dimen¬ 
sions and at right boundary where the flow is supersonic all primitive variables p, u, v and p are extrapolated with first 
order approximation. Eigure 11 shows the pressure contours using 40 x 40 Q4 mesh. 


8.4.2. Oblique Shock Reflection from aflat plate: 

In this test case 137] the domain is rectangular [0, 3] x [0, 1]. The boundary conditions are 


1. Inflow (left boundary ): p - I, u - 2.9, v — 0,p- 1/1.4 

2. Post shock condition (top boundary) : p - 1.69997, u - 2.61934, v = -0.50633, p = 1.52819 

3. Bottom boundary is a solid wall where slip boundary condition is applied, i.e., \.n - 0. 

4. At right boundary where the flow is supersonic all primitive variables p, u, v and p are extrapolated with first 
order approximation. 


Pressure plots for 60 x 20, 120 x 40 and 240 x 80 quadrilateral mesh are given in figure [T^ The comparison of residue 
plots are given in figure 13 Eor triangular unstructured mesh (number of nodes: 2437 and number of triangles: 4680) 
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Figure 6: Strong Rarefactions Riemann Problem 




Figure 7: Solution of 2D Burgers equation on Q4 mesh 


the pressure contours are given in figure 14 and residue plot is shown in figure 15 
are captured quite accurately at correct positions. 


The incident and reflected shocks 


8.4.3. Supersonic flow over a half cylinder: 

Two supersonic test cases with inflow Mach numbers 2 and 3 are tested on a half cylinder |[39l. The domain is 
half circular, the left outer circle is inflow boundary. Small circle inside the domain is a cylinder wall and the straight 
edges on right sides are supersonic outflow boundaries. Pressue plots (see figure [T6] l show that the bow shock in front 
of the half-cylinder is captured accurately at the right position in each case which are compared with existing results 


20 



























X X 

Figure 8: Solution of 2D Burgers equation on T3 mesh 
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No. of Iterations 


Figure 9: Residue plot for 2D Burgers equation on Q4 mesh 



X 



Figure 10: Solution of Burgers equation with original SUPG method using Delta parameter (Left) and the proposed 
method without any Delta parameter (Right) on 32 x 32 Q4 Mesh. 


ia. 

S.5. Numerical Experiments for Implicit KSUPG 

All previously solved ID test cases are again solved for implicit KSUPG method. 
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Figure 11: Pressure contours (0.18:0.005:0.29) using 40 x 40 Q4 elements. 


8.5.1. Sod’s Shock Tube Problem: 

The number of node points are 100 and CFL number is 0.6. Final time is t = 0.01. Figure 
velocity, pressure and mach number plots. 


17 shows the density, 


8.5.2. Shock Tube Problem of Lax: 

The number of node points are 100 and CFL number is 0.6. Final time is t=0.13. Figure 
velocity, pressure and internal energy plots. 


18 shows the density. 


8.5.3. Strong Rarefactions Riemann Problem: 

The number of node points are 200 and CFL number is 0.6. Final time is t=0.15. Figure 
pressure and velocity plots. 


19 shows the density. 


8.6. Comparison of Explicit and Implicit KSUPG scheme for 2D Euler test case 

In case of implicit KSUPG scheme, 2D Euler test cases are solved using CFL =1. For the comparison of explicit 
and implicit KSUPG schemes shock reflection test case is solved using different grids. Figure shows the pressure 
contour plot on 120 x 40 Q4 grid and the residue vs number of iterations. 

Iteration speed-up ratio is defined as the ratio of number of iterations required for the residue to drop below 
a predefined tolerance value for an explicit method to that of an implicit method. Similarly, one can define the 
computational speed-up ratio which is the ratio of total computational time required for explicit method to that of 
implicit method. Following tables shows comparison of computational cost and number of iterations taken for explicit 
and implicit KSUPG schemes for oblique shock reflection test case. 



Total Computational Cost 

Grid Size 

60x20 

120x40 

240 X 80 

Explicit 

KSUPG 

8388 sec 

69948 sec 

189000 sec 

Implicit 

KSUPG 

1338 sec 

8784 sec 

37288 sec 

Computational 
Speed-up Ratio 

6.26 

7.96 

5.06 
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Figure 13; Residue plots for 60 x 20, 120 x 40 and 240 x 80 



Figure 14: Pressure contours (0.8:0.1:2.8) using T3 element. 


8.6.1. Oblique shock 

Oblique shock test case is solved using implicit KSUPG method. Figure shows the pressure contours using 
40 X 40 Q4 mesh. 

8.7. Sparsity Pattern in Explicit and Implicit KSUPG Method 

Figure shows the sparsity patterns for explicit and implicit KSUPG method for shock reflection test case over 
60 X 20 mesh size using Q4 elements. Both matrices are unsymmetric. The number of nonzero entries in explicit 
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Figure 15; Residue plots for a triangular elements. 



X X 


Figure 16: Pressure contours (0.6:0.2:3.6) for Mach 2 (left) and (1:0.5:7) for Mach 3 (right) using 46x46 quadrilateral 
mesh. 

method is 41296 and in case of implicit method it is 145100 which are far more than the previous case. The half band¬ 
width of explicit KSUPG matrix is 245 whereas for implicit KSUPG it is 248. In case of implicit KUPG method the 
condition number of coefficient matrix of assembled system is high compared to explicit KSUPG scheme. Following 
table shows the condition number calculated using L 2 norm for different grid size for shock reflection test case. 
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Figure 17: Sod’s Shock Tube Problem 



Condition Number 

Grid Size 

Implicit KSUPG (CFL =1) 

Explicit KSUPG 

60x20 

2.1670e04 

8.3722e03 

120 X 40 

9.4669e04 

3.4646e04 

240 X 80 

4.18822e05 

1.40931e05 


9. Conclusions 

In this paper we presented a novel explicit as well as implicit kinetic theory based streamline upwind Petrov 
Galerkin scheme (KSUPG scheme) in finite element framework for both scalar case (inviscid Burgers equation in ID 
and 2D) and vector case (ID and 2D Euler equations of gas dynamics). The proposed numerical scheme is simple 
and easy to implement. The important advantage in using a kinetic scheme in finite element framework is that, 
instead of dealing with nonlinear hyperbolic conservation laws, one needs to deal with a simple linear convection 
equation. In comparison with the standard SUPG scheme, the advantage of the proposed scheme is that, it does not 
contain any complicated expression for t (which is the intrinsic time scale) especially for vector equations. Also, for 
the multidimensional Burgers equation, standard SUPG scheme requires additional diffusion term (shock capturing 
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Figure 18: Shock Tube Problem of Lax 


parameter) which is not required in the proposed scheme. The accuracy and robustness of the scheme is demonstrated 
by solving various test cases for Burgers equation and Euler equations. Spectral stability analysis is done for 2D 
linear equation which gives an implicit expression of stable time step. Finally, comparison between explicit and 
implicit versions of KSUPG scheme is done with respect to the number of iterations, computational cost, sparsity 
pattern and the condition number of a global system of equations. 
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